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Abstract: The article derives multivariate Generalized Gram-Charlier (GGC) 
series that expands an unknown joint probability density function (pdf) 
of a random vector in terms of the differentiations of the joint pdf of a 
reference random vector. Conventionally, the higher order differentiations 
of a multivariate pdf in GGC series will require multi-element array or 
tensor representations. But, the current article derives the GGC series in 
vector notations. The required higher order differentiations of a multivari¬ 
ate pdf in vector notations are achieved through application of a specific 
Kronecker product based differentiation operator. Overall, the article uses 
only elementary calculus of several variables; instead Tensor calculus; to 
achieve the extension of an existing specific derivation for GGC series in uni¬ 
variate to multivariate. The derived multivariate GGC expression is more 
elementary as using vector notations compare to the coordinatewise ten¬ 
sor notations and more comprehensive as apparently more nearer to its 
counterpart for univariate. The same advantages are shared by the other 
expressions obtained in the article; such as the mutual relations between 
cumulants and moments of a random vector, integral form of a multivariate 
pdf integral form of the multivariate Hermite polynomials, the multivariate 
Gram-Charlier A (GCA) series and others. 
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1. Introduction 

The Generalized Gram-Charlier (GGC) series expands an unknown pdf as a 
linear combination of the increasing order differentiations of a reference pdf, 
where the coefficients of expansion involve cumulant differences between those 
of an unknown pdf and a reference pdf. The GGC expansions are used to ap¬ 
proximate pdf and functions of pdf in Statistics, Machine Learning, Economics, 
Chemistry, Astronomy and other application areas. There have been used Pois¬ 
son’s distribution (Aroian, 1937), log-normal distribution (Aroian, 1947), bino¬ 
mial distribution (Rietz, 1927), gamma distribution (Bowers, 1966) and oth¬ 
ers as reference pdfs. But, the Gram-Charlier (GC) expansion with Gaussian 
density as a reference pdf is the most popular and identified as the Gram- 
Charlier A (GCA) series. Specifically, the GCA series is used for test of Gaus- 
sianity or near Gaussian pdf approximations (Aroian, 1947; Viehland, 1994; 
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Girolami and Fyfe, 1996; Jondeau and Rockinger, 2001; Del Brio et al., 2009), 
for entropy measure and independence measure approximations (Amari et al., 
1996; Hyvarinen et al., 2001), for optima analysis through derivatives of pdf 
(Boscolo and Roychowdhury, 2003), for time-frequency analysis (Cohen, 1998) 
and others. The rearrangement of the terms in GCA series results into the 
Edgeworth series with better convergence property. 

There exists multiple ways to derive univariate GCA series, as reported by 
Hald (2000); Hald and Steffensen (2002). It’s generalization to univariate GGC 
series is derived by Schleher (1977); Cohen (1998); Berber an-Santos (2007); 
Cohen (2011). 

The multivariate GGC or GCA series derivation requires multivariate 
representations of the Taylor series, the increasing order differentiations of a 
reference pdf and the cumulants. To signify the representation issue in required 
multivariate extensions; it is worth quoting Terdik (2002) that says, ‘ ... though 
the generalizations to higher dimensions may be considered straightforward, 
the methodology and the mathematical notations get quite cumbersome when 
dealing with the higher order derivatives of the characteristic function or other 
functions of a random vector ... ’. 

Conventionally, the higher order differentiations of a multivariate pdf; and 
therefore, the multivariate cumulants and multivariate Hermite polynomials; 
require multilinear representations. It is known and acknowledged historically 
that going from matrix like notations to tensor notations for multivariate cumu¬ 
lants and multivariate Hermite polynomials have made the representation more 
transparent and the calculations more simpler (McCullagh, 1987). Though the 
tensor notations have advantages over the matrix notations, they require com¬ 
ponentwise separate representations and are more tedious compare to the vector 
notations. As a unihed and comprehensive solution to this, there has been used 
an approach based only on elementary calculus of several variables by Terdik 
(2002). The approach uses a specific Kronecker product based differential opera¬ 
tor, identified as the ‘K-derivative operator’, to achieve vectorization of the Jaco¬ 
bian matrix of a multivariate vector function. The successive applications of this 
K-derivative operator achieves vectorization of the higher order derivatives also. 
Using this approach, there have been derived the multivariate Taylor series and 
the higher order cumulants (Terdik, 2002; S. Rao Jammalamadaka and Terdik, 
May, 2006); as well, the multivariate Hermite polynomials (Terdik, 2002) in vec¬ 
tor notations resulting into more transparent representations. In fact, it could 
be noticed that the same approach has been first used to derive vector Her¬ 
mite polynomials by Holmquist (1996); and later^ it has been formalized and 
generalized by Terdik (2002). 

There exists various approaches deriving multivariate GCA series with 
various representations. They include GCA series representation using multi¬ 
element matrix notations for moments and cumulants by Sauer and Heydt (1979); 
using recursive formula for Hermite polynomials by Berkowitz and Garner (1970); 
using tensor representation for cumulants and Hermite polynomials by McCullagh 

^Somehow, the citation for (Holmquist, 1996) is not found in article (Terdik, 2002). 
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(1987, Chapter 5); using tensor representation for cumulants and involving mul¬ 
tivariate Bell polynomials by Withers and Nadarajah (2014); using vector mo¬ 
ments and vector Hermite polynomials by Holmquist (1996) and others. There 
exists various derivations for multivariate Edgeworth series (Davis, Dec., 1976; 
Skovgaard, 1986; Amari and Kumon, 1983; McCullagh, 1987; Kollo and von Rosen, 
1998; Withers and Nadarajah, 2014). There also exists multivariate GGC series, 
derived by McCullagh (1987, Chapter 5), in tensor notations. But, as per the 
author’s knowledge, there exists neither the multivariate GGC series nor the 
multivariate Edgeworth series in vector notations. Eor the ease of the readers 
in following and comparing the various representations; the existing represen¬ 
tations of multivariate GGC series and multivariate GCA series are shown and 
compared in A. 

Overall, to take advantages due to the recent advancement in represen¬ 
tation, this article extends a specific derivation for univariate GGC series by 
Berberan-Santos (2007) to multivariate; using only elementary calculus of sev¬ 
eral variables instead of Tensor calculus. As a by product, it also derives mutual 
relations between vector cumulants and vector moments of a random vector; 
integral form of the multivariate pdf; integral form of the multivariate vector 
Hermite polynomials and the multivariate GCA series. All the derived multi¬ 
variate expressions are more elementary as using vector notations and more 
comprehensive as apparently more nearer to their counterparts for univariate; 
compare to their coordinatewise tensor notations. The intermediate theoretical 
results, in the article, are verified using suitable known examples. 

Towards the aim of the article, the next Section 2 briefs some neces¬ 
sary background on the Kronecker product and a way to obtain vectorization 
of the higher order differentiations of a multivariate pdf. It also obtains the re¬ 
quired multivariate Taylor series expansion using the derived notations. After 
the preliminary background, this article follows almost the same sequence for 
multivariate as that in (Berberan-Santos, 2007) for univariate. The Section 3 
uses the characteristic function and the generating functions to derive cumulants 
and moments of a random vector in vector notations with their mutual relation¬ 
ships. The Section 4 obtains multivariate pdf in terms of its vector cumulants. 
The expressions for derivatives of multivariate Gaussian density and vector Her¬ 
mite polynomials are derived in Section 5. The Section 6 derives multivariate 
GCA series by representing an unknown pdf in terms of the Gaussian pdf as 
a reference. The Section 7 derives GGC Expansion, representing an unknown 
pdf in terms of the a known reference pdf The Section 8 derives an unknown 
characteristic function of a random vector in terms of a reference characteristic 
function. The Section 9 derives the same GGC expansion in a more compact 
way that summarizes the approach of the whole derivation. Finally, Section 10 
concludes the article. For the sake of clarity; the calculation details, the proofs 
and the expressions for existing multivariate expansions are kept in appendix 
at the end of the article. 
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2. Vectorization of the higher order differentiations 

The section briefs the Kronecker product and the way it can be applied to 
achieve vectorization of the higher order differentiations of a multivariate pdf. 
Based on it, the multivariate Taylor series is obtained in vector notations. More 
details can be found on the Kronecker Product in (Magnus and Nendecker, 
1999, Chapter 2), on achieving vectorization of the higher order differentia¬ 
tions in (Terdik, 2002; S. Rao Jammalamadaka and Terdik, May, 2006) and on 
the commutation matrices in (Magnus and Nendecker, 1999, Chapter 3, Section 

7 ). 

Definition 1 (Kronecker Product Operator (®)). The Kronecker Product Op¬ 
erator (®) between matrices A with size pxq and B with size mxn is defined as: 


aiiB 

ai2B • 

■ aigB 

a2iB 

a22B 

■ a2gB 

OpiB 

ap2B 

UpgB 


( 1 ) 


The resultant matrix is of dimension pm x qn. As a further example; let A 
is with size p x 1 and B is with size m x 1, then A 0 B' is^ a matrix with 
size p X m. A 0 A is symbolically represented as A®^ and has size p^ x In 
general, A 0 A 0 ... 0 A (n times) is symbolically represented as A®" and has 
size p"' X 1. 

Definition 2 (Jacobian Matrix). Let A = (Ai, A 2 ,..., \d)', A € and f(A) = 
(/i(A), / 2 (A), ..., fmWY G R™ be CL differentiable m-component vector function. 
Then, Jacobian matrix o/f(A) (J(f)) is an m x d matrix defined as under: 




■ af 5f di ' 

d\i ’ 8X2 ’ " ' ’ dXd 


dfi dh 
9Ai 9 A 2 

dh 

9Ai 


dim dU 

9 Ai 9A2 


dh 1 
d\d 


( 2 ) 


dfm 
d\d . 


Let the vector differential operator be defined as a column vector Dx = 




/ 9 9 9 

V9Ai ’ 9A2 ’ • ■ • ’ dXd ) 

re-written as: 


then the Jacobian matrix, in terms of the D\, can be 


J(f(A)) = DA(f) = f(A)D'x = (/i(A), / 2 (A),..., fm{\))' 


/ d d 



This implies that to match the definition of differentiation from matrix calculus, 
the vector differential operator should be applied from the right to the left. This 

^The symbol ’ stands for Transpose of a matrix 
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is same as the requirement to be satisfied on generalization of vector derivative to 
matrix derivative as discussed by Magnus (2010). So, applying vector derivative 
operator from right to the left, has been kept as a rule throughout the article. 


Definition 3 (The K-derivative Operator). Let A = (Ai, A 2 ,..., A^)', A € 
the vector differential operator Dx = gfj) ■ • ■ i 917 ) ® differentiable 

m-component vector function f(A) = (/i(A), / 2 (A),..., fmW)' G Then, the 
Kronecker product between Dx and f(A) is given as under: 







r ^ 

dh . 

. Ml. 1 


■ /i(A) ■ 


1- -1 

O 

d\\ 


dXi 

9A2 

9 Ad 

Dff(A) = 

/ 2 (A) 

0 

d 

0X2 

= Vec 

df 2 

dXi 




. /m(A) . 


d 

- dXd - 


d.fm 

dX\ 

dU . 
9A2 

. Mn. 
dXd . 


=.D®f(A) = rec(§)' = v-,,(|.f) 


(4) 

(5) 


where, the V ec operator converts mx d matrix into an md x 1 column vector by 
stacking the columns one after an other. The operator D® is called Kronecker 
derivative operator or simply, K-derivative operator. 


Thus, the Kronecker product with the vector differential operator, obtains 
vectorization of the transposed Jacobian of a vector function. Corresponding to 
the definition, the order differentiation is given by: 

Df f = Df (Df-If) = [/i(A), / 2 (A),..., fm{X)]' 0 

( 6 ) 


d d d 

dXi ’ 9A2’ ’ dXd 


The is a column vector of dimension mdf x 1. Some important properties 

of the K-derivative operator, those are useful in the further derivations, are 
listed in Appendix B. 


2.1. Application of the K-derivative operator to the multivariate 
Taylor series 

Let X = {Xi, X 2 ,..., Xd)' be a d-dimensional column vector and /(x) be the 
function of several variables differentiable in each variable. Using the defined 
K-derivative operator, the Taylor series for /(x), expanding it at origin, is given 
as: 

m—oo ^ 

/(x)= ^ —c(m,d)'x®- (7) 

ml 

m—O 

where, c(m, d) is the vector of dimension d™ x 1 and given in terms of the 
derivative vector Dx = afl’ ‘ ‘ ‘ ’ sfj) 

c(m,d) = (Dr/(x))|x=o (8) 
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The Taylor series expansion near A = 0, called the Maclaurian series, of some 
required functions based on the Equation (7) are derived in appendix C. 

3. Moments, cumulants and characteristic function of a random 
vector 

Let X = (All, ^ 2 ,be a d-dimensional random vector and /(x) be its 
joint pdf differentiable in each variable. 

The Characteristic function (J^) of x is defined as the expected value of 
eix where A = (Ai, A 2 ,..., Ad)^ A G Also, both the characteristic function 
and the pdf are the Fourier Transform (F) of each other, in the sense they are 
dual. 

J-x(A) = A{e*(-'^)| = F(/(x)) = r ... r /(x)e*("'^)dx (9) 

^ J —00 J — 00 

Expanding e™ using its Maclaurian series in Equation (80) in appendix C, we 
get: 


-7^x(A) = 


( 10 ) 


where, m(fc, d) is the order moment vector of dimension x 1 and given by 

m(fe,d)= [ x®V(x)dx 
Js.<‘ 

( 11 ) 


Also, /(x) = F ^(J'(A)) = F M ^m(fc,d) 

\fe =0 

m(fc, d)' / 1 


k\ 


= E' 


00 




k\ (27r)^ 

^g^_l),m(^DWd(x) 

fc =0 

(■.■ Proof in Appendix D.l) 

The Moment Generating Function (MCE) of /(X) is given as 
M(A) = = / /(X)e’''^dX 




( 12 ) 


= ^m(fc,d) 


,A® 


fc =0 


fc! 


( •.■ Expanding e’^ 


Assuming M(A) and J^{X) are expanded using Taylor series, 

ni{k,d) = Df M(A)|x=o = (-i)'=Df Jx(A)|a=o 


(13) 

(14) 

(15) 
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The Cumulant Generating Function (CGF) of /(X) is given by, 

C(A) = lnM(A)==^c(fc,d)'— (16) 

k=l 

where, c(fc, d) is the order cumulant vector of dimension x 1. 

The Cumulant Generating Function (CGF) of /(X) can also be defined using 
the characteristic function, as under: 


C(A)=lnJ-(A)=^c(fc,d)'li^ (17) 

/t -1 

Assnming C(A) and C(A) have been expanded using Taylor series, 

c(fc, d) = Df'=C(A)|A^o = (-*)"Df'=C x(A)|a=o (18) 

Significantly, this section has derived the moments and the cumulants of a ran¬ 
dom vector in vector notations. 


3.1. Relation between the cumulant vectors and the moment vectors 


The relation between the moments and the cumulants is given by combining 
Equation (14) and Equation (16) as below: 

oo -(gifc / oo 

M(A) = y] m(fe, d)'— = exp I y] c(fe, d)'— \ (19) 

k—0 \fc—1 / 


For fc = 1, using Equation (15), we get: 

m(l,d) =Df M(A)k=o 


Applying K-derivative (D®) to Equation (19), 


m(p, d) 

p^l 




(p- 1)! 


A=0 = y^c(p,d) 


^ 0 l/ 


) exp 


Y^(d,dy 


(p-l)! 

A®«' 


\g=l 


<?! 


A=0 


m(l, d) = c(l, d) 


Similarly, based on Equation (15), for any k taking order K-derivative of 
above Equation (19) on both sides relates m(fc,d) with c{k,d). For example, 
the cases for fe = 2 and fe = 3 are shown in Appendix D.2. Overall, the vector 
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moments in terms of the vector cumulants can be summarized as under: 


m(l, d) = c(l, d) 
m(2, d) = c(2, d) + c(l, 

m(3, d) = c(3, d) + 3c(2, d) 0 c(l, d) + c(l, d)®^ 
m(4, d) = c(4, d) + 4c(3, d) 0 c(l, d) + 3c(2, d)®^ + 6c(2, d) 
0 c(l,d)®2+c(l,d)®4 

m(5, d) = c(5, d) + 5c(4, d) 0 c(l, d) + 10c(3, d) 0 c(2, d) 

+ 10c(3, d) 0 c(l, d)®^ + 15c(2, d)®^ 0 c(l, d) 

+ lOc(2,d)0c(l,d)®3 + c(l,d)®5 
m(6, d) = c(6, d) + 6c(5, d) 0 c(l, d) + 15c(4, d) 0 c(2, d) 

+ 15c(4, d) 0 c(l, d)®^ + 10c(3, d)®^ + 60c(3, d) 0 c(2, d) 
0 c(l, d) + 20c(3, d) 0 c(l, d)®2 + 15c(2, d)®^ 

+ 45c(2, d)®2 0 c(l, d)®2 + 15c(2, d) 0 c(l, d)®^ + c(l, d)®® 


The above set of equations can be represented through more compact formulas as 


under. The Equation (21) gives generalized order d-variate cumulant vector 


in terms of the moment vectors and the Equation (22) gives the vice-a-versa. 



( 21 ) 


P—0 

1 



( 22 ) 


where, is a specific commutation matrix with corresponding dimensions 


that changes the place of the cumulants for Kronecker product such that the 


expression has decreasing order cumulants from left to the right, i.e. the higher 
order cumulant vector on the left and the lower order cumulant vector on the 
right. As the Kronecker products are non-commutative, without using the com¬ 
mutation matrices it would have been impossible to derive the compact formula. 
The derived multivariate expressions reduce to the following expressions for di¬ 
mension d = 1 and are exactly same as those derived in (Berberan-Santos, 
2007). 



(23) 



Thus, the derived multivariate expressions in Equation (20) are elementary vec¬ 
tor extensions to those for univariate. 
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4. Multivariate pdf representation in terms of the cumulants 


From Equation (11) and Equation (17), the multivariate pdf /(x) can be written 
as: 


/(X) = F-‘(e««) = exp(-ix'A)dA 

(24) 

As pdf is a real function and Re(e"^+®'®e“*‘") = cos{B — C) = e^ cos(C' — B), 
the Equation (24) can be re-written as: 


/(x) = 


(27r)“ 


exp 




\k^l 


cos l^x'A + £ £^-^(*)2fc(A)®(2'=-i)^ dA 
The integrand in this equation is an even function. So, 

J(R+r 2fc! j 

, y e(2>:-l d)' ^ 


COS X 


(25) 


(26) 


where, R’*' = {x € R : a: > 0}. The Equation (25) and the Equation (26) give 
a multivariate pdf in terms of the cumulants. As they are derived using Taylor 
series expansion, the infinite order differentiability is an implicit assumption. 

The equations can be verified using known pdf examples with finite num¬ 
ber of moments and cumulants. Let say, the impulse delta density function has 
only the first order cumulant being non-zero and all other higher order cumu¬ 
lants are zero. Using this knowledge in Equation (25), 

/ cos((x-c(l,d))'A) dA (27) 

JR<* 

= d(x — c(l, d)) (shifted impulse delta function) 


Let’s take another example, the Gaussian density function has first two order 
cumulants nonzero and all other order cumulants are zero. Using this knowledge 
in Equation (25), 

= jir L ((=4 - <=(1. -i))' " ,28, 

= G(x) ( See Appendix D.3 for proof) 
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5. The multivariate Hermite polynomials in integral form 

An interesting application of the integral form of multivariate pdf representation 
is achieved in this section. The multivariate Gaussian expressed as in Equation 
(28) is used to derive it’s differentials and Hermite polynomials in a simple way. 
Taking K-derivative of G(x), 

Gik) Df'=G(x) = exp 

cos ^(x — c(l, d))^ A + c?A (29) 

The multivariate Hermite polynomials defined by Holmquist (1996) are defined 
as under: 


Hfc(x; 0, Cx) = [G(x; 0, Cx)]-\-l)'= (CxD,)®'= G(x; 0, Cx) (30) 
where, G(x;0,Cx) = ICxp^^^ (27r)“‘’*/^ exp ^-ix'C,^^x^ 

= ICxp'/" (2^)-‘’*/2 exp {Vec C"')' x®^^ (31) 

This is equivalent to the 1-dimensional definition of Hermite polynomials^ by 
Rodrigues’s formula in Equation (32), except the introduction of matrix Cx. 

Hk{x) = [Gix)]~^{-l)'^-^Gix) where, G{x) = (32) 


Using Equation (29), the Equation (30) for multivariate Hermite polynomials is 
rewritten as: 



Taking Cx = Id, where Id is dxd identity matrix; using the property {Vec Id)'x®^ 
x'x and using change of variable as \l\/2 = u - the integral form of multivariate 
Hermite polynomials is obtained as under: 


Hfe(x;0,ld) = (2) 2 (tt) ‘’'G exp (-x'x 


f u®^ exp (— u'u) cos f v^x'u — du (34) 

\ 2 / 


3 This is the ’probabilists’ Hermite polynomials and not the ’physicists’ Hermite polyno¬ 
mials used by Berberan-Santos (2007). 
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The result in Equation (29) can also be obtained using Equation (28) and 
applying the derivative property of Fourier transform (F). The derivative of 
G(x) is given by, 

(x) = F-i (f(G('=) (x))) (35) 

= 79^/ (*'’^)®''F(G(x))exp-*’‘'^dA (36) 

jR’i 

' L cos (^(x - c(l.d))' A + f) dA 

(37) 


6. Multivariate Gram-Charlier A series 


Till now, the article has derived - an unknown pdf expressed in terms of its 
cumulants in Equation (25); the Gaussian density function expressed in terms 
of its cumulants in Equation (28) and the Hermite polynomials in Equation 
(34). Based on them, the multivariate Gram Gharlier A series that expresses 
an unknown pdf using Gaussian density as a reference can be obtained. The 
expansion assumes first and second order cumulants being same for both the 
unknown pdf and the reference pdf. Using the expansion exp(A+i?) cos(G+I?) = 
exp(A) exp(i?)(cos G cos D — sin G sin D), the Equation (25) can be re-written 
as: 


|cos((x - c(l, d))'\) cos £^^^(*)2'=(A)^(2 '=-i)^ 
- sin((x - c(l, d))'\) sin 


d\ (38) 


Using the expansions in Appendix G, parts of the Equation (38) can be approx¬ 
imated upto maximum 6*^-order statistics as under: 


exp 




2k\ ^ (2 fc-l)! 

cfe =2 / \ fc =2 ^ ' 


= 1-F 


c(4,d)'A®4 c(6,d)'A®®\ 1 /c(3,d)'A®3\®^ 


4! 


6 ! 


3! 


(39) 


exp 




c(3,d)'A'^^ c(5,d)'A 


/ \ (S>5 


3! 


5! 


(40) 
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Using above Equation (39) and Equation (40), the Equation (38) can be re¬ 
written as: 

cos((x-c(l,d))'A) 

(l + - i (c(6, d) - 10c(3, A®= + ...) I dX 

■ (2^ L cos ((X - c(l,<i))'A - 

c(3^_c(5^^ (41, 


Using the Equation (29) for derivatives of Gaussian defined, the above Equation 
can be simplified as: 


/(x) = G(x) - h|^G<»>(x) + d(l2lG<-‘>(x) - d|<G<“>(x) 

_^ c(6,d)' + 10c(3,d)»^ GW(x) + ... (42) 

6 ! 

The Equation (42) is the Gram-Charlier A series expressed directly in terms 
of the cumulants and the derivatives of the Gaussian pdf. Usually, the GGA 
is represented in terms of the Hermite polynomials. So, the GGA expansion 
(Equation (42)) in terms of the Hermite polynomials; either using definition in 
Equation (30) or using Hfc(x; 0, C~^) derived in Equation (33); can be re-written 
as: 


/(x) = G(x) 


c(3, d)' 


1 + (C-')®' H3(x; 0, C,,) + (C;')”" H4(x; 0, C^) 


c(4,d)' 


3! 


+ £AA (cp)®' H 5 (x; 0, CJ + <=(<i.'lY + 10c(3,dr' 


/(x) = G(x) 


5! 

H 6 (x; 0 , Cx) -l- ... 

c(3, d)' 


1 -h 


3! 


(43) 

H3(x; 0, Id) + £i4^H4(x; 0, Id) + £%^H5(x; 0, Id) 


c(6, d)' + 10c(3, d)®^ 


6 ! 


4! 


H6(x;0,Id) -f ... 


5! 


(44) 


Einally, the GGA series, in vector notations, can be expressed either using 
order derivative of Gaussian (Gfc(x)) or using order vector Hermite polyno¬ 
mials (Hfc(x)) as under: 

/(x)=f;(-l)'=^i^G«(x) (45) 

= f:^Hdx;0,l4) (46) 
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7. Multivariate Generalized Gram-Charlier series 


To derive the generalized Gram-Charlier series, an unknown pdf /(x) need be 
represented in terms of any known reference pdf tfix), where both the pdfs 
are represented in terms of their cumulants. Let the order cumulant vector 
of the reference pdf be Cr{k,d). Then, the order cumulant difference 
vector S{k, d) is: 6{k, d) = c(fc, d) — Cr{k, d),yk. Using 6{k, d), the Equation (25) 
can be re-written as under: 

= (2^ L (I (g 




COS I x'A -t I 


'^fc=l 


Cr{2k — 1, d) 

{2k-l)\ 




(47) 


12^ L (£ (f 


cos X 


'A-f ^ 


Cr{2k - l,dy 2k 


- sin g'A + ^ 


d(2fc- l,d)' 

(2fc- 1)! 

fe-i) 1 


(48) 
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Using the expansions in Appendix C, parts of the Equation (48) can be approx¬ 
imated upto maximum 6 *^-order statistics as under: 


exp 




\k^l 


\k^l 


[2k-l)\ 


=1+1: 


5{2k, d)'!. ^ ',02fc 


Kk=l 


2k\ 


{i\f 




5(2k, d)' ^ ^(g,2fc 


02 






2k\ 

02 


(*A)' 


= 1 + - 


S{2,d)'^^ , .5(4, d)' 


2 ! 


4! 


^04_ d(6, d) 

6 ! 


1 / d{ 2 ,dy . 


2 \ 2 ! 

03 




02 


d)' > 0)2 

2 ! 


d(4, d)' |g |4 
4! 


(d(l,d)'A)®^ 


d-Y >02 

2 ! 


02 


0(d(l,d)'A)®" 


+ ^ A®^) - 1 (-d(l, dyxf^ - 0 d(l, d)'A®4^ +... 

= 1 - i (d(l, d)®" + d( 2 , d))' A ®2 + 1 (d(i, d)®4 5^(2, d) 0 d(l, d )®2 

+3d(2, d)®" + 4d(3, d) 0 d(l, d) + d(4, d))' A®4 - ^ (d(l, d)®® + 15d(2, d)®® 

+ 10d(3, d)®2 + 15d(4, d) 0 d(2, d) + 15d(4, d) 0 d(l, d)®^ -h 20d(3, d) 0 d(l, d)®^ 
+15d(2, d) 0 d(l, d)®4 + 45d(2, d)®^ 0 d(l, d)®^ 6d(5, d) 0 d(l, d) 

+60d(3, d) 0 d( 2 , d) 0 d(l, d) + d( 6 , d))' A®® -h ... (49) 
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Similarly, 


exp 




= (-^( 1 , dy\ + 


^(2,d)'A®2 5(l,d)'A 


2 ! 


^(3,d)'A®3 ^ ^(2,d)'A®2 ^ ^(4,d)'A®‘‘ ^ 5(l,d)'A®i 


3! 2! 4! 

1 /^( 2 ,d)'A®^\®^ ^(l,d)'A®^ 


2 ! 


1 ! 


- (-<5(l,d)'A)®^ + ... 


1 ! 


= ^( 1 , dyx + i (5(3, d) + 35(2, d) 0 5(1, d) + 5(1, d)®3)' a®3 _ 1 (^( 5 , d) 
+55(4, d) 0 5(1, d) + 155(2, d)®^ ^ ^( 1 ^ d) + 105(3, d) 0 5(2, d) 

+105(2, d) 0 5(1, d)®3 + 105(3, d) 0 5(1, d)®^ + 5(1, d)®*^)' A®® ... (50) 

Now, D®^'0(x) can be obtained by taking fc*^-order K-derivative of Equation 
(25) as under; 

^(fe)(x) = DfV(x) = j^Jj\r'^exp c,(2M)' (^A)^2fc^ 


COS 


(“x'A + D + bli JA (51) 


The above Equation (51) with the previous results on expansions in Equation 
(49) and Equation (50) can be used to simplify the Equation (48). This derives 
the Generalized Gram-Charlier (GGC) series expressing an unknown pdf /(x) in 
terms of the cumulant difference vectors (5(fc, d)) and derivatives of a reference 
pdf tp^^yx) as under: 


/(x) = f;(-l)^^(|^V'('=)(x) (52) 

/c =0 
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where, 


a(0, d) = 1 

a(l, d) = ^(1, d) 

a(2,d) = ^(2,d) + ^(l,d)®2 

a(3, d) = 6(3, d) + 36(2, d) ® ^(1, d) + ^(1, 

a(4, d) = d(4, d) + 4d(3, d) 0 d(l, d) + 36(2, d)®^ + 6d(2, d) 

0(5(l,d)®2 + d(l,d)®^ 

q:( 5, d) = d(5, d) + 5d(4, d) 0 d(l, d) + 10d(3, d) 0 <5(2, d) (53) 

+ 10d(3, d) 0 <5(1, d)®2 + 15^(2, d)®2 0 <5(1, d) 

+ 10d(2, d) 0 d(l, d)®3 + d(l, d)®5 
a(6, d) = d(6, d) + 6d(5, d) 0 d(l, d) + 15d(4, d) 0 <5(2, d) 

+ 15d(4, d) 0 d(l, d)®2 + iod(3, d)®2 + 60<5(3, d) 0 d(2, d) 

0 d(l, d) + 20d(3, d) 0 <5(1, d)®2 + 15<5(2, d)®^ + 45^(2, d)®^ 

0 d(l, d)®2 + 15<5(2, d) 0 d(l, d)®4 + d(l, d)®® 

The above set of equations (53) has exact resemblance with that expressing mo¬ 
ments in terms of the cumulants in Section 3.1. This must happen, as Equation 
(52) for GGC expansion with <5(x) as a reference pdf is matching Equation (12). 
This matching proves that ck(k, d) is related in same way to 6(k, d), as m(fc, d) 
to c(k, d). That is,: 



(54) 


Further, the a(k,d) in Equation (52) recursively can be obtained in terms of 
the cumulant difference vector 6(k,d) as under: 



(55) 


p—0 


where, is a specific commutation matrix; as described previously; to 

change the order of the cumulants for Kronecker product such that the expres¬ 
sion has decreasing order cumulants from left to the right. 

The verification of the derived GGC can be obtained by taking Gaussian 
density as a reference pdf. With Gaussian density as a reference, d(l,d) = 0 
and 6(2, d) = 0. So, the coefficients a(k,d) in Equation (52) can be derived as 
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under: 

q:( 0, d) = 1 
q:(1 , d) = 0 
q:( 2, d) = 0 

cy.(3,d) = d{3,d) = c{3,d) (56) 

q:(4, d) = ^(4, d) = c(4, d) 

q:(5, d) = ^(5, d) = c(5, d) 

a(6, d) = S(6, d) + 10^(3, = c(6, d) + 10c(3, 

Thus, the GGC series is derived and verified using known examples. 


8. Characteristic function of an unknown random vector in terms of 
a reference characteristic function 

The GGC derived as in Equation (52) can be used to give the characteristic 
function of an unknown pdf, in terms of the characteristic function of a reference 
pdf. For that taking Fourier transform (F) of Equation (52), we get: 


J-x(A) = f;(-l)^-^%^F(V^«(x)) 


k—0 


^x(A)= 


^Q:(fc,d) 


fc! 

,(*A)« 


(57) 


. fc =0 


k\ 


F(^(x)) (■.■ differentiation property of F ) 

(58) 


= exp 


^d(fe,d) 


,(*A)« 


.k=l 


k\ 


Tripd) (■.• Equation (54)) 


(59) 


where, is the characteristic function of the reference pdf 


9. Compact derivation for the Generalized Gram-Charlier expansion 

The compact derivation of Equation (52) follows as under: 

,,(a)®'=' 


-Fx(A) = exp 


^c(fc,d) 


= exp 


oo 


Y,s{k,d) 


.k=l 


kl 

,(fA)®'= 

k\ 


( definition in Equation (17) ) (60) 

(iA)®'= 


exp 


^Crik,d)- 


.k=l 


kl 


('.■ S{k, d) = c{k, d) — Cr{k, d)) 

,{t\r 

a) 

. fc =0 


^a{k,dy 


kl 


F(V’(x)) 


(61) 


(62) 
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Taking inverse Fourier transform of the above equation brings 

/_ 1 \k 

/(x) = * V'(x) (63) 

/t-0 

or/(x) = J V’^'=)(x) (64) 

k—0 

where, * indicates convolution. Thus, the Equation (52) is obtained in a more 
compact way. 

10. Conclusion 

The article has derived multivariate Generalized Gram-Charlier (GGG) expan¬ 
sion in Equation (52); combined with Equation (55); that expresses an unknown 
multivariate pdf in terms of vector cumulants and vector derivatives of a refer¬ 
ence pdf. The multivariate Gram-Gharlier A series is derived in Equation (45) 
and Equation (46) representing an unknown multivariate pdf in terms of its 
vector cumulants and vector Hermite polynomials. There has been also derived 
compact formulas for obtaining multivariate vector moments from vector cu¬ 
mulants in Equation (21) and vise-a-verse in Equation (22); the integral form 
of multivariate pdf representation in Equation (25) and the integral form of 
multivariate vector Hermite polynomials in Equation (33), as well, in Equation 
(34). The expressions are derived using only elementary calculus of several vari¬ 
ables in vector notations through Kronecker product based derivative operator. 
Thus, they are more transparent and more comprehensive compare to their 
corresponding multi-linear matrix representations or tensor representations. 


Appendix A: The Multivariate Representations of GCA Series and 
GGC Series 


As mention in Section 1, this section of the appendix describes existing repre¬ 
sentations of GCA series and GGC series. The goal is to place together various 
historical representations for the ease of comparison, on the level of difficulty or 
simplicity in representation, to the readers. Therefore, no attempt is made to 
explain their derivation or each terms in representation. For further details the 
actual references need be referred. 

The GCA series representation using multi-element matrix notations for 
cumulants and moments by Sauer and Heydt (1979) is as under: 


OO OO OO 


d 


/x(x) = ^ ^ • 

E 

Cs,s,-sv ■ (-1)^--*^ 

(65) 

Si—0 S2—0 

Sd=0 

L J 


with, CsiS 2 . 

■Sd — 


(66) 
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where, is the constant depending upon cross-moments and Hi(x) is 

the one-dimensional Hermite polynomial of order. 

The GCA Series representation using recursive formula for Hermite poly¬ 
nomials by Berkowitz and Garner (1970) is as under: 


/x(x) = G(x) ^ 
m—0 


with, Am 



Jm{z)f:^{x)dx 


(67) 

( 68 ) 


where, G(x) denote the multivariate Gaussian; z = (x — /.i) is the normalized 
variable; /x is the mean vector; {i7m(x)} and {Jm(x)} are the complete bi- 
orthogonal system of Hermite polynomials. 

Using recursive relations of {Jm(x)}, Am is given recursively as under: 


Am. — 


N 




. 2=1 


N 

2=1 


N 

2=1 


A 

■^m —Cfc —e 


,k = 1, . . . ,d 


(69) 


where, N is the number of available samples and e*, is a vector with a ”1” as 
the component and ” 0” elsewhere. This defines the coefficients of expansions 
also recursively. 

The GGC series representation using tensor notations for cumulants and 
Hermite polynomials by McCullagh (1987, Chapter 5) is as under: 

/x(a:;K) = /o(x) [l -k 77*/ii(x) -k (x)/2! -k 77*'^'"/iyfc(x)/3! -b ...] (70) 


where, hi{x) = /i,(x; A) = /i(x)//o(x), hij{x) = hij{x-X) = /y(x)//o(x), ... 
and /i(x) = dfo{x)/dx\fij{x) = d'^foix.)/dx ^,...; so on. Also, given 
K^, ■ ■ ■ are the cumulant tensors of random vector x and A*, ... 

are the cumulant tensors of the reference pdf /o(x); we get: 


2 2 \ 2 2,2 i.j \i 

T] = K — X = K — A 




- A^ 




The formal ‘moments’ rf ^ ■ ■ ■ s-re defined based on the formal ‘cumu¬ 
lants’ (or the cumulant differences) rf gQ 

Taking /o(x) = G(x) i.e. multivariate Gaussian density as the reference 
pdf and taking 77* = 0 ,r]^d = q in above Equation (70); the GCA series based 
on cumulant tensors is written as under: 


/x(x; k) = G(x) [1 + s:*’^'’'=/i,,fe(x)/3! + 

+ + 10 ^^,Xk^l,m,nj . . .] (71) 

As could be observed, the GGC series and GCA series using tensor notations 
adds quite an ease to representation. But, with increase in number of terms, the 
difficulty in representation increases. 
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The GCA series using vector moments and vector Hermite polynomials 
by Holmquist (1996) is as under: 


/x(x) = G(x - ^Ji■ Cx) ^ -Cl (x - /.; C^)E {Hfe(x - /x; C^} (72) 

fc =0 

where, Gfc(x—/x) is the order vector derivative of G(x—/x) and E {Hfe(x — /x)} 
is the expectation of order vector Hermite polynomial that is the function 
of vector moments. 


Ef{nk{x 


/x;C-i)®^}=fc!S 


dl. 


[fc/2] 

E 

j=0 


Tna-k- 2 j G {-Vec Cx)®-^ 
{k-2j)lj& 


(73) 


where, £1/ {(X — /x)®-^ } = mj(/x). 

The GCA series representation using Bell polynomials is obtained by 
Withers and Nadarajah (2014). Here, the Bell polynomials are represented through 
cumulant tensors. With r = (ri, r 2 ,..., rd) € I'*, x*" = r! = ri!r 2 ! ■ ■ .rdl 

and |r| = n + r 2 + ... + r^; the GCA series is as under: 

OO 

/x(x)/G(x) - 1 = ^ B,i£,(x, Cx)/r! (74) 

|r|>3 

|r| 

where, Brj (75) 

3=0 

B-'.x = E Q (j' + 1 - kl+lBr-l, (76) 

kr = K{Xrl, Xr2, ■ ■ ■ , Xrd) (77) 

where, kr is the order cumulant and k ( x ) is the moment. 

Finally, the existing representations described in this Section A can be 
compared with those derived in the current article. , can be compared with 
the GCA series derived in Equation (45) and the GGC series derived in Equa¬ 
tion (52) combined with Equation (54). More specifically, the GCA series in 
Equation (65) using multi-element array representation and the GGC series in 
Equation (70) using tensor representation can be compared with the GCA series 
derived in Equation (45) and the GGC series derived in Equation (52), com¬ 
bined with Equation (54), in vector notations. The comparison demonstrates 
the ease obtained in representation. The same advantages are also obtained for 
other intermediate results in the article. 


Appendix B: Some important properties of the K-derivative 
operator 

Some important properties of the K-derivative operator are listed below. 
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Property 1 (Scaling Property). Let A = (Ai, A 2 ,..., A^)', A S f(A) £ 
and fi(A) = Af(A), where A is an n x m matrix. Then 

Df(fi) = (A®Irf)Df(f) (78) 

where, 1 ^ is a d-dimensional unit matrix. 

Property 2 (Chain Rule). Let A = (Ai, A 2 ,..., A^)', A £ f(A) £ and 
g(A) £ R™C Then 

Df(f 0 g) = d)((Df f) 0 g) + f 0 (Df g) (79) 

where, K 3 <_). 2 (mi, m 2 , d) is a commutation matrix of size mim 2 d x mim 2 d that 
changes the order of the the Kronecker product components. For example, 

K3^2(mi,m2, d)(ai 0 a2 0 as) = ai 0 as 0 a2 

The above K-derivative properties can be used to derive further the following 
properties: 

1. Df A = Pec Id 

2. Df A®'' = k 0 Vec Id) 

3. Let g{X) = sl' f{X) 

Df 5 (A) = (a 0 ld)'Df/(A) 

4. Df a'A®'' = fc (a 0 Id)' 0 Vec Id) = k a' 0 Id) 

5. Df a'®'=A®'= = k (a®'= 0 Id)' (^A®'^"^ 0 Vec Id) = k (a'A)®''"^ 0 a 

6 . Df’'A®'= = fc(fc - 1 ) • • • (fc - r + 1 ) (^A®''-’’ 0 {Vec Id)®’') 

( Repeated application of the Chain Rule ) 

7. Df’'a'®'=A®'' = k{k-l) ■ ■ ■ (fc-r+1) (a®'= 0 I®’')' (^A®'="’' 0 {Vec Id)®’’) 

= k{k — 1) • • • (fc — r + 1) (a'A)®^”’’ 0 a®*" 

( Repeated application of the Chain Rule and the property: ) 
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Appendix C: The Taylor series expansion of some required 
functions near zero 


The Taylor series expansion of the required functions near A = 0, based on the 
Equation (7), are given as under: 


,,, ^ (a'x)®'= , (a'x)®2 (a'x)®3 

= - —TT— = 1 + a'x + -h ^ ^ 


k=0 

OO 


k\ 


2 ! 


3! 


sill' 


fe =0 


■ (I \ (-1)^ ( , \mk+i / (a'x)'^^ (a'x)®^ 

in(a'x) = > -^(a'x)'^"^'"+^ = a'x - ^f-h -r— 

^ ^ ^ 2fc + i r 3 5 




(80) 


(81) 

(82) 


e“'^ cos(b'y) = 1 + a'x + 


(a'x )®2 (b'y )®2 (a'x)®3 (a'x) 0 (b'y)®^ 


2! 2! 3! 

(a'x)®4 (b'y)®4 (a'x )®2 ^ (b'y)®^ 


2 ! 


4! 


e"^ sin(b'y) = b'y + (a'x) 0 (b'y) — 


4! 2!2! 

(a'x )®2 0 (b'y) (b'y)®3 


(83) 


2! 3! 

(a'x) 0 (b'y)®3 (b'y)®5 (a'x)®^ ^ (b'y)®^ 


3! 


5! 


2!3! 


+ ... (84) 


Appendix D: Some Proofs 
D.l. K-derivative of d{x.) 

Based on the Differentiation Property of Fourier Transform, the following is 
obtained: 

F-i (F(D®V(x)) = F-i ((*A)®''F(/(x))) 

D®'=5(x) = F-i(a)®'= 
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D.2. Relationship between m(fc, d) and c{k, d) for k — 2 and k — 3 


For fe = 2; using Equation (15), Equation (19) and the differentiation rules for 
Kronecker product in Appendix B; the following is derived: 


i(2,d) = DfM(A)|A=o 


p=2 


{p-2)\ 


(85) 

A-O = I 




®(p-2) 


' Id ® Id 


(p-2)! 


/ oo A®?' 


A=0 


+ “i 1X] 

I P=1 




®2 


(p- 1)! 

^m(2, d) = c(2, d) + c(l, d)®^ 


/ °° \®g' 


A=0 


( 86 ) 


where, p3 O 2 € tp is the required permutation and Kp 3 <_). 2 (di, ^ 2 , d) is the 
corresponding commutation matrix. 

For fc = 3; using Equation (15), Equation (19) and the differentiation rules for 
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Kronecker product in Appendix B; the following is derived: 

m(3,d) =DfM(A)|A=o (87) 


p=3 




(p- 3)! 


=0 — Kp3^2(^l) ^ 2 , d)Kp3^2(^li ^^ 2 , d) 




- 3 ) 0 1^3 


,p=3 


ip-3)1 


'exp '^c{q,dy 


+ K:p3V^2(di,d2,d) ^ |^c(p,d) 

i p=2 


\q=l 




A=0 


(p-2)! 


\p=i 


(a®(p-i) 0 Id) ' 


(p- 1)! 


/ 00 x«ig' 


v=o + Kp3'^2 (di, d2, d) 


2 I J2 

I p=2 


^(A®(^’-2)0ld2)' 


®1 




(p-2)! 


^c(p,d) 


^ (a®(p-i) 0 Id) ' 


(g)l 


\ 


(p- 1)! 


/ 00 x®g ' 

>exp ( ^c{q,dy — 


\q=l 


01. 


A=0 


^c{p, d) 

I P=1 


^ (a®(p-i) 0 Id) 


02 


(p- 1)! 


(gii 


/ 00 x®g' 


iP-^)- ^ 

=4>m(3, d) = c(3, d) + 3c(2, d) 0 c(l, d) + c(l, d)®^ 


A=0 


^p=l 


( 88 ) 


D.3. Multivariate Gaussian representation in terms of the 
cumulants (The proof related to the Section 4 

Applying D® on Equation (28), 

Df/(x) = (A 0 Id) 0 exp sin ((x - c(l, d))' A) dA 


(89) 

Let, U{\) = exp ^- L.I_J_a®^^ (90) 

D®;7(A) = -c(2, d)' (A 0 Id) exp (91) 

Also, let V (A) = (x — c(l, d)) cos ((x — c(l, d))' A) (92) 

^ f EdA = sin ((x — c(l, d))^ A) (93) 
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Using above Equation (90), Equation (92) in Equation (89) and taking inv c(2, d) 
Vec we get: 

Dx /(x) = £ (Df C/(A)) (^J Vd/j dX (94) 

Applying integration by parts to above Equation (94), we get: 

Df/(x) = |^-ic(2,d)'A®^^ sin ((x - c(l, d))'A) 

— exp (x — c(l, d)) cos ((x — c(l, d))' A) dA| (95) 

= inn c(2, d)'(x — c(l, d))/(x) (96) 

The solution of above differential equation leads to the following: 

/(x) = cexp (— {inv c(2, d))' (x — c(l, d))) (97) 

where, c = /(O) = exp 

= (27r)-''/2 (|C,|)-V2 (98) 

^ /(x) = (27r)“'^/^ (|Cx|)“^^^ exp (- {inv c(2, d))' (x - c(l, d))) (99) 

^ (2^ Jd A®^^ cos (x'A - c(l, d)'A) dA = G(x) (100) 
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